package scu.maqiang.fes;

import static scu.maqiang.numeric.Constants.GP3;
import static scu.maqiang.numeric.Constants.GW3;

import java.util.Arrays;

import scu.maqiang.mesh.Mesh;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.CG;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.Pair;
import scu.maqiang.numeric.ParamCheck;
import scu.maqiang.numeric.SRMatrix;

public class FES3H201 extends FES {

	public FES3H201(Mesh mesh) {
		super(mesh, 1);
		eleStiff = FES3H201::elementStiff;
		eleMass = FES3H201::elementMass;
		eleSource = FES3H201::elementSource;
	}

	public double getValue(double[] u, double[] xyz) {
		ParamCheck.checkEqual(u.length, this.ndof);
		int[] eleIdx = new int[this.nPerEle];
		double[] eleVal = new double[this.nPerEle];
		return getValue(u, xyz, eleIdx, eleVal);
	}

	public double getValue(double[] u, double[] xyz, int[] eleIdx, double[] eleVal) {
		ParamCheck.checkEqual(u.length, this.ndof);
		ParamCheck.checkEqual(eleIdx.length, this.nPerEle);
		ParamCheck.checkEqual(eleVal.length, this.nPerEle);
		double[] uu = new double[20];
		double value = 0;
		double[][] coord = new double[20][3];
		for (int i = 0, n = mesh.getNt(); i < n; i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coord);
			if (CG.pointInH8(coord, xyz)) {
				CG.H20ShapeFunction(coord, xyz, eleVal);
				for (int j = 0; j < 20; j++) {
					uu[j] = u[ele[j]];
					eleIdx[j] = ele[j];
					value += uu[j] * eleVal[j];
				}
				return value;
			}
		}
		throw new IllegalArgumentException("The Point ( " + Arrays.toString(xyz) + " ) is not in the domain!");
	}

	public Pair<Double, Integer> getValue(double[] u, double[] xyz, int previousEleIdx, double[][] coef, int[] eleIdx,
			double[] eleVal) {
		ParamCheck.checkEqual(u.length, this.ndof);
		ParamCheck.checkEqual(eleIdx.length, this.nPerEle);
		ParamCheck.checkEqual(eleVal.length, this.nPerEle);
		double value = 0;
		double x = xyz[0], y = xyz[1], z = xyz[2];
		double[] v = { 1.0, x, y, z, x * x, y * y, z * z, x * y, y * z, z * x, x * x * y, x * y * y, y * y * z,
				z * z * y, x * x * z, z * z * x, x * y * z, x * x * y * z, x * y * y * z, x * y * z * z };
		double[][] coord = new double[20][3];
		int[] ele = mesh.getElement(previousEleIdx);
		mesh.getCoordPerElement(ele, coord);
		if (CG.pointInH8(coord, xyz)) {
			// CG.H20ShapeFunction(coord, xyz, eleVal);
			MVO.matmul(coef, v, eleVal);
			for (int j = 0; j < 20; j++) {
				// eleIdx[j] = ele.idx[j];
				value += u[ele[j]] * eleVal[j];
			}
			// System.out.println("Yes!");
			return new Pair<Double, Integer>(value, previousEleIdx);
		} else {
			for (int i = 0, n = mesh.getNt(); i < n; i++) {
				ele = mesh.getElement(i);
				mesh.getCoordPerElement(ele, coord);
				if (CG.pointInH8(coord, xyz)) {
					// CG.H20ShapeFunction(coord, xyz, eleVal);
					CG.H20ShapeFunctionCoef(coord, coef);
					MVO.matmul(coef, v, eleVal);
					for (int j = 0; j < 20; j++) {
						eleIdx[j] = ele[j];
						value += u[ele[j]] * eleVal[j];
					}
					return new Pair<Double, Integer>(value, i);
				}
			}
			throw new IllegalArgumentException("The Point ( " + Arrays.toString(xyz) + " ) is not in the domain!");
		}

	}

	/**
	 * �õ�����Ԫ�����ض�������µ�ֵ
	 * 
	 * @param u      ����Ԫ�����飬����Ϊ����Ԫ�����ɶ���
	 * @param xyz    ������������
	 * @param data   Pair���󣬵�һ������Ϊ����Ԫ���ڸ������µĺ���ֵ���ڶ�������õ������������ڵĵ�Ԫ���
	 * @param coef   �õ�
	 * @param eleIdx
	 * @param eleVal
	 */
	public void getValue(double[] u, double[] xyz, Pair<Double, Integer> data, double[][] coef, int[] eleIdx,
			double[] eleVal) {
		ParamCheck.checkEqual(u.length, this.ndof);
		ParamCheck.checkEqual(eleIdx.length, this.nPerEle);
		ParamCheck.checkEqual(eleVal.length, this.nPerEle);
		double value = 0;
		double x = xyz[0], y = xyz[1], z = xyz[2];
		double[] v = { 1.0, x, y, z, x * x, y * y, z * z, x * y, y * z, z * x, x * x * y, x * y * y, y * y * z,
				z * z * y, x * x * z, z * z * x, x * y * z, x * x * y * z, x * y * y * z, x * y * z * z };
		double[][] coord = new double[20][3];
		int[] ele = mesh.getElement(data.getSecond());
		mesh.getCoordPerElement(ele, coord);
		if (CG.pointInH8(coord, xyz)) {
			// CG.H20ShapeFunction(coord, xyz, eleVal);
			MVO.matmul(coef, v, eleVal);
			for (int j = 0; j < 20; j++) {
				// eleIdx[j] = ele.idx[j];
				value += u[ele[j]] * eleVal[j];
			}
			// System.out.println("Yes!");
			data.setFirst(value);
			return;
			// return new Pair<Double, Integer>(value, previousEleIdx);
		} else {
			for (int i = 0, n = mesh.getNt(); i < n; i++) {
				ele = mesh.getElement(i);
				mesh.getCoordPerElement(ele, coord);
				if (CG.pointInH8(coord, xyz)) {
					// CG.H20ShapeFunction(coord, xyz, eleVal);
					CG.H20ShapeFunctionCoef(coord, coef);
					MVO.matmul(coef, v, eleVal);
					for (int j = 0; j < 20; j++) {
						eleIdx[j] = ele[j];
						value += u[ele[j]] * eleVal[j];
					}
					data.setFirst(value);
					data.setSecond(i);
					return;
				}
			}
			throw new IllegalArgumentException("The Point ( " + Arrays.toString(xyz) + " ) is not in the domain!");
		}

	}

	public void getValue_Cube(double[] u, double[] xyz, Pair<Double, Integer> data, double[][] coef, int[] eleIdx,
			double[] eleVal) {
		ParamCheck.checkEqual(u.length, this.ndof);
		ParamCheck.checkEqual(eleIdx.length, this.nPerEle);
		ParamCheck.checkEqual(eleVal.length, this.nPerEle);
		double value = 0;
		double x = xyz[0], y = xyz[1], z = xyz[2];
		double[] v = { 1.0, x, y, z, x * x, y * y, z * z, x * y, y * z, z * x, x * x * y, x * y * y, y * y * z,
				z * z * y, x * x * z, z * z * x, x * y * z, x * x * y * z, x * y * y * z, x * y * z * z };
		double[][] coord = new double[20][3];
		int[] ele = mesh.getElement(data.getSecond());
		mesh.getCoordPerElement(ele, coord);
		if (CG.pointInCube(coord, xyz)) {
			// CG.H20ShapeFunction(coord, xyz, eleVal);
			MVO.matmul(coef, v, eleVal);
			for (int j = 0; j < 20; j++) {
				// eleIdx[j] = ele.idx[j];
				value += u[ele[j]] * eleVal[j];
			}
			// System.out.println("Yes!");
			data.setFirst(value);
			return;
			// return new Pair<Double, Integer>(value, previousEleIdx);
		} else {
			for (int i = 0, n = mesh.getNt(); i < n; i++) {
				ele = mesh.getElement(i);
				mesh.getCoordPerElement(ele, coord);
				if (CG.pointInCube(coord, xyz)) {
					// CG.H20ShapeFunction(coord, xyz, eleVal);
					CG.H20ShapeFunctionCoef(coord, coef);
					MVO.matmul(coef, v, eleVal);
					for (int j = 0; j < 20; j++) {
						eleIdx[j] = ele[j];
						value += u[ele[j]] * eleVal[j];
					}
					data.setFirst(value);
					data.setSecond(i);
					return;
				}
			}
			throw new IllegalArgumentException("The Point ( " + Arrays.toString(xyz) + " ) is not in the domain!");
		}

	}

	/**
	 * ���㵥Ԫ�նȾ���
	 *
	 * @param ele
	 * @param xyz
	 * @param kxx
	 * @param tp
	 * @param eleMatrix
	 */
	public static void elementStiff(int[] ele, double[][] xyz, double[][] kxx, BVPType tp, double[][] eleMatrix) {
		double[] x = MVO.col(xyz, 0);
		double[] y = MVO.col(xyz, 1);
		double[] z = MVO.col(xyz, 2);
		double[] dx = new double[20];
		double[] dy = new double[20];
		double[] dz = new double[20];
		double[] v = new double[20];
		double xi, et, zt, detJ, wxarea, c;
		MVO.fill(eleMatrix, 0.0);
		switch (tp) {
		case COMMON:
			for (int i = 0; i < 3; i++) {
				xi = GP3[i];
				for (int j = 0; j < 3; j++) {
					et = GP3[j];
					for (int k = 0; k < 3; k++) {
						zt = GP3[k];
						detJ = CG.ISOMap3D(CG::H20ShapeFunction, x, y, z, xi, et, zt, v, dx, dy, dz);
						c = MVO.dot_product(kxx[0], v);
						wxarea = GW3[i] * GW3[j] * GW3[k] * detJ * c;
						for (int si = 0; si < 20; si++) {
							for (int sj = 0; sj < 20; sj++) {
								eleMatrix[si][sj] += (dx[si] * dx[sj] + dy[si] * dy[sj] + dz[si] * dz[sj]) * wxarea;
							}
						}
					}
				}
			}
			break;
		default:
			throw new IllegalArgumentException("Wrong BVPType!");
		}

	}

	public static void elementMass(int[] ele, double[][] xyz, double[][] coef, BVPType tp, double[][] eleMatrix) {
		double[] x = MVO.col(xyz, 0);
		double[] y = MVO.col(xyz, 1);
		double[] z = MVO.col(xyz, 2);
		double[] dx = new double[20];
		double[] dy = new double[20];
		double[] dz = new double[20];
		double[] v = new double[20];
		double xi, et, zt, detJ, wxarea, c;
		MVO.fill(eleMatrix, 0.0);
		switch (tp) {
		case CONSIST_MASS:
			for (int i = 0; i < 3; i++) {
				xi = GP3[i];
				for (int j = 0; j < 3; j++) {
					et = GP3[j];
					for (int k = 0; k < 3; k++) {
						zt = GP3[k];
						detJ = CG.ISOMap3D(CG::H20ShapeFunction, x, y, z, xi, et, zt, v, dx, dy, dz);
						c = MVO.dot_product(coef[0], v);
						wxarea = GW3[i] * GW3[j] * GW3[k] * detJ * c;
						for (int si = 0; si < 20; si++) {
							for (int sj = 0; sj < 20; sj++) {
								eleMatrix[si][sj] += v[si] * v[sj] * wxarea;
							}
						}
					}

				}
			}
			break;
		case LUMP_MASS:
			for (int i = 0; i < 3; i++) {
				xi = GP3[i];
				for (int j = 0; j < 3; j++) {
					et = GP3[j];
					for (int k = 0; k < 3; k++) {
						zt = GP3[k];
						detJ = CG.ISOMap3D(CG::H20ShapeFunction, x, y, z, xi, et, zt, v, dx, dy, dz);
						c = MVO.dot_product(coef[0], v);
						wxarea = GW3[i] * GW3[j] * GW3[k] * detJ * c;
						for (int si = 0; si < 20; si++) {
							eleMatrix[si][si] += v[si] * wxarea;
						}
					}

				}
			}
			break;
		default:
			throw new IllegalArgumentException("Wrong BVPType!");
		}
	}
	
	public void assembleHeatThetaScheme(ScalarFunc[] coefFunc, double[] param, BVPType[] tp, SRMatrix LMatrix, SRMatrix RMatrix) {
		double[][] kxx = new double[1][nPerEle];
		double[][] rho = new double[1][nPerEle];
		double[][] eleMatrixStiff = new double[dofPerElement][dofPerElement];
		double[][] eleMatrixMass = new double[dofPerElement][dofPerElement];
		double[][] eleLMatrix = new double[dofPerElement][dofPerElement];
		double[][] eleRMatrix = new double[dofPerElement][dofPerElement];
		double dt = param[0];
		double theta = param[1];
		for (int i = 0; i < nE; i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coordRegion);
			int label = mesh.getElementLabel(i);
			for (int j = 0, nj = nPerEle; j < nj; j++) {
				kxx[0][j] = coefFunc[0].action(coordRegion[j], label, param);
				rho[0][j] = coefFunc[1].action(coordRegion[j], label, param);
			}
			elementStiff(ele, coordRegion, kxx, tp[0], eleMatrixStiff);
			elementMass(ele, coordRegion, rho, tp[1], eleMatrixMass);
			for(int ii = 0; ii < dofPerElement; ii++ ) {
				for(int jj = 0; jj < dofPerElement; jj++) {
					eleLMatrix[ii][jj] = eleMatrixMass[ii][jj] + theta * dt * eleMatrixStiff[ii][jj];
					eleRMatrix[ii][jj] = eleMatrixMass[ii][jj] - (1 - theta) * dt * eleMatrixStiff[ii][jj];
				}
			}
			LMatrix.assemble(ele, eleLMatrix);
			RMatrix.assemble(ele, eleRMatrix);
		}
	}
	
	public void assembleHeatThetaScheme(double[] coefFunc, BVPType[] tp, SRMatrix LMatrix, SRMatrix RMatrix) {
		double[][] kxx = new double[1][nPerEle];
		double[][] rho = new double[1][nPerEle];
		double[][] eleMatrixStiff = new double[dofPerElement][dofPerElement];
		double[][] eleMatrixMass = new double[dofPerElement][dofPerElement];
		double[][] eleLMatrix = new double[dofPerElement][dofPerElement];
		double[][] eleRMatrix = new double[dofPerElement][dofPerElement];
		Arrays.fill(kxx[0], coefFunc[0]);
		Arrays.fill(rho[0], coefFunc[1]);
		double dt = coefFunc[2];
		double theta = coefFunc[3];		
		for (int i = 0; i < nE; i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coordRegion);
			elementStiff(ele, coordRegion, kxx, tp[0], eleMatrixStiff);
			elementMass(ele, coordRegion, rho, tp[1], eleMatrixMass);
//			System.out.println(MVO.toString(eleMatrixMass));
//			System.out.println(MVO.toString(eleMatrixStiff));
//			System.exit(0);
			for(int ii = 0; ii < dofPerElement; ii++ ) {
				for(int jj = 0; jj < dofPerElement; jj++) {
					eleLMatrix[ii][jj] = eleMatrixMass[ii][jj] + theta * dt * eleMatrixStiff[ii][jj];
					eleRMatrix[ii][jj] = eleMatrixMass[ii][jj] - (1 - theta) * dt * eleMatrixStiff[ii][jj];
				}
			}
			LMatrix.assemble(ele, eleLMatrix);
			RMatrix.assemble(ele, eleRMatrix);
		}
	}

	public static void elementSource(int[] ele, double[][] xyz, double[][] f, BVPType tp, double[] eleF) {
		double[] x = MVO.col(xyz, 0);
		double[] y = MVO.col(xyz, 1);
		double[] z = MVO.col(xyz, 2);
		double[] dx = new double[20];
		double[] dy = new double[20];
		double[] dz = new double[20];
		double[] v = new double[20];
		double xi, et, zt, detJ, wxarea, c;
		MVO.setValue(eleF, 0.0);
		switch (tp) {
		case COMMON:
			for (int i = 0, n = GP3.length; i < n; i++) {
				xi = GP3[i];
				for (int j = 0; j < n; j++) {
					et = GP3[j];
					for (int k = 0; k < n; k++) {
						zt = GP3[k];
						detJ = CG.ISOMap3D(CG::H20ShapeFunction, x, y, z, xi, et, zt, v, dx, dy, dz);
						c = MVO.dot_product(f[0], v);
						wxarea = GW3[i] * GW3[j] * GW3[k] * detJ * c;
						for (int si = 0; si < 20; si++) {
							eleF[si] += v[si] * wxarea;
						}
					}
				}
			}
			break;
		default:
			throw new IllegalArgumentException("Wrong BVPType!");
		}
		// System.out.println(Arrays.toString(ele));
	}

	public void computeGradient2(double[] u, double[][] gradU) {
		int len = u.length;
		ParamCheck.checkEqual(len, this.ndof);
		ParamCheck.checkEqual(len, gradU[0].length);
		ParamCheck.checkEqual(len, gradU[1].length);

		double volume;
		double[][] coord = new double[20][3];
		double[][] grad = new double[20][3];
		double[] volumeCount = new double[this.GetNdof()];
		double[] eleU = new double[20];
		for (int i = 0, n = mesh.getNt(); i < n; i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coord);
			MVO.getSlice(u, ele, eleU);
			CG.gradientValue3DH20(coord, eleU, grad);
			volume = CG.volume3DH8(coord);
			for (int j = 0; j < 20; j++) {
				gradU[0][ele[j]] += grad[j][0] * volume;
				gradU[1][ele[j]] += grad[j][1] * volume;
				gradU[2][ele[j]] += grad[j][2] * volume;
				volumeCount[ele[j]] += volume;
			}
		}
		for (int i = 0, n = mesh.getNv(); i < n; i++) {
			gradU[0][i] /= volumeCount[i];
			gradU[1][i] /= volumeCount[i];
			gradU[2][i] /= volumeCount[i];
		}
	}

	public void computeSecondDerivative(double[] u, double[][] sdU) {
		int len = u.length;
		ParamCheck.checkEqual(len, this.ndof);
		ParamCheck.checkEqual(len, sdU[0].length);
		ParamCheck.checkEqual(len, sdU[1].length);
		ParamCheck.checkEqual(len, sdU[2].length);

		double volume;
		double[][] coord = new double[20][3];
		double[][] eleSDU = new double[20][6];
		double[] angleCount = new double[this.GetNdof()];
		double[] eleU = new double[20];
		for (int i = 0, n = mesh.getNt(); i < n; i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coord);
			for (int j = 0; j < 20; j++) {
				eleU[j] = u[ele[j]];
			}
			CG.secondDerivativeValue3DH20(coord, eleU, eleSDU);
			volume = CG.volume3DH8(coord);
			for (int j = 0; j < 20; j++) {
				sdU[0][ele[j]] += eleSDU[j][0] * volume;
				sdU[1][ele[j]] += eleSDU[j][1] * volume;
				sdU[2][ele[j]] += eleSDU[j][2] * volume;
				sdU[3][ele[j]] += eleSDU[j][3] * volume;
				sdU[4][ele[j]] += eleSDU[j][4] * volume;
				sdU[5][ele[j]] += eleSDU[j][5] * volume;
				angleCount[ele[j]] += volume;
			}
		}
		for (int i = 0, n = mesh.getNv(); i < n; i++) {
			sdU[0][i] /= angleCount[i];
			sdU[1][i] /= angleCount[i];
			sdU[2][i] /= angleCount[i];
			sdU[3][i] /= angleCount[i];
			sdU[4][i] /= angleCount[i];
			sdU[5][i] /= angleCount[i];
		}

//		for(int i = 0; i < mesh.getNt(); i++) {
//			ele = mesh.getElement(i);
//			for(int k = 0; k < 6; k++) {
//				for(int j = 0; j < 20; j++) {
//					eleU[j] = sdU[k][ele.idx[j]];
//				}
//				sdU[k][ele.idx[8]] = 0.5 * (eleU[0] + eleU[1]);
//				sdU[k][ele.idx[9]] = 0.5 * (eleU[1] + eleU[2]);
//				sdU[k][ele.idx[10]] = 0.5 * (eleU[2] + eleU[3]);
//				sdU[k][ele.idx[11]] = 0.5 * (eleU[3] + eleU[0]);
//				sdU[k][ele.idx[12]] = 0.5 * (eleU[4] + eleU[5]);
//				sdU[k][ele.idx[13]] = 0.5 * (eleU[5] + eleU[6]);
//				sdU[k][ele.idx[14]] = 0.5 * (eleU[6] + eleU[7]);
//				sdU[k][ele.idx[15]] = 0.5 * (eleU[7] + eleU[4]);
//				sdU[k][ele.idx[16]] = 0.5 * (eleU[0] + eleU[4]);
//				sdU[k][ele.idx[17]] = 0.5 * (eleU[1] + eleU[5]);
//				sdU[k][ele.idx[18]] = 0.5 * (eleU[2] + eleU[6]);
//				sdU[k][ele.idx[19]] = 0.5 * (eleU[3] + eleU[7]);
//			}
//		}
	}

	public void computeFirstAndSecondDerivative(double[] u, double[][] fdU, double[][] sdU) {
		int len = u.length;
		ParamCheck.checkEqual(len, this.ndof);
		ParamCheck.checkEqual(len, sdU[0].length);
		ParamCheck.checkEqual(len, sdU[1].length);
		ParamCheck.checkEqual(len, sdU[2].length);

		double volume;
		double[][] coord = new double[20][3];
		double[][] eleFDU = new double[20][3];
		double[][] eleSDU = new double[20][6];
		double[] angleCount = new double[this.GetNdof()];
		double[] eleU = new double[20];
		for (int i = 0, n = mesh.getNt(); i < n; i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coord);
			for (int j = 0; j < 20; j++) {
				eleU[j] = u[ele[j]];
			}
			CG.firstAndSecondDerivativeValue3DH20(coord, eleU, eleFDU, eleSDU);
			volume = CG.volume3DH8(coord);
			for (int j = 0; j < 20; j++) {
				fdU[0][ele[j]] += eleFDU[j][0] * volume;
				fdU[1][ele[j]] += eleFDU[j][1] * volume;
				fdU[2][ele[j]] += eleFDU[j][2] * volume;
				sdU[0][ele[j]] += eleSDU[j][0] * volume;
				sdU[1][ele[j]] += eleSDU[j][1] * volume;
				sdU[2][ele[j]] += eleSDU[j][2] * volume;
				sdU[3][ele[j]] += eleSDU[j][3] * volume;
				sdU[4][ele[j]] += eleSDU[j][4] * volume;
				sdU[5][ele[j]] += eleSDU[j][5] * volume;
				angleCount[ele[j]] += volume;
			}
		}
		for (int i = 0, n = mesh.getNv(); i < n; i++) {
			fdU[0][i] /= angleCount[i];
			fdU[1][i] /= angleCount[i];
			fdU[2][i] /= angleCount[i];
			sdU[0][i] /= angleCount[i];
			sdU[1][i] /= angleCount[i];
			sdU[2][i] /= angleCount[i];
			sdU[3][i] /= angleCount[i];
			sdU[4][i] /= angleCount[i];
			sdU[5][i] /= angleCount[i];
		}
	}

}
